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In a first step toward the comprehension of neural activity, one should focus on the 
stability of the possible dynamical states. Even the characterization of an idealized regime, 
such as that of a perfectly periodic spiking activity, reveals unexpected difficulties. In this 
paper we discuss a general approach to linear stability of pulse-coupled neural networks 
for generic phase-response curves and post-synaptic response functions. In particular, we 
present: (1) a mean-field approach developed under the hypothesis of an infinite network 
and small synaptic conductances; (2) a "microscopic" approach which applies to finite but 
large networks. As a result, we find that there exist two classes of perturbations: those 
which are perfectly described by the mean-field approach and those which are subject to 
finite-size corrections, irrespective of the network size. The analysis of perfectly regular, 
asynchronous, states reveals that their stability depends crucially on the smoothness 
of both the phase-response curve and the transmitted post-synaptic pulse. Numerical 
simulations suggest that this scenario extends to systems that are not covered by the 
perturbative approach. Altogether, we have described a series of tools for the stability 
analysis of various dynamical regimes of generic pulse-coupled oscillators, going beyond 
those that are currently invoked in the literature. 



Keywords: linear stability analysis, splay states, synchronization, neural networks, pulse coupled neurons, Floquet 
spectrum 



1. INTRODUCTION 

Networks of oscillators play an important role in both bio- 
logical (neural systems, circadian rhythms, population dynam- 
ics) (Pikovsky et al, 2003) and physical contexts (power grids, 
Josephson junctions, cold atoms) (Hadley and Beasley, 1987; 
Filatrella et al., 2008; Javaloyes et al, 2008). It is therefore com- 
prehensible that many studies have been and are still devoted to 
understanding their dynamical properties. Since the development 
of sufficiently powerful tools and the resulting discovery of gen- 
eral laws is an utterly difficult task, it is convenient to start from 
simple setups. 

The first issue to consider is the model structure of the single 
oscillators. Since phases are typically more sensitive than ampli- 
tudes to mutual coupling, they are likely to provide the most 
relevant contribution to the collective evolution (Pikovsky et al., 
2003). Accordingly, here we restrict our analysis to oscillators 
characterized by a single, phase-like, variable. This is typically 
done by reducing the neuronal dynamics to the evolution of the 
membrane potential and introducing the corresponding veloc- 
ity field which describes the single-neuron activity. Equivalently, 
one can map the membrane potential onto a phase variable and 
simultaneously introduce a phase-response curve (PRC) [Upon 
changing variables, the velocity field can be made independent of 
the local variable (as intuitively expected for a true phase). When 
this is done, the phase dependence of the velocity field is moved 
to the coupling function, i.e., to the PRC] to take into account 
the dependence of the neuronal response on the current value of 
the membrane potential (i.e., the phase). In this paper we adopt 



the first point of view, with a few exceptions, when the second one 
is mathematically more convenient. 

As for the coupling, two mechanisms are typically invoked 
in the literature, diffusive and pulse-mediated. While the former 
mechanism is pretty well understood [see e.g., the very many 
papers devoted to Kuramoto-like models (Acebron et al., 2005)], 
the latter one, more appropriate in neural dynamics, involves a 
series of subtleties that have not yet been fully appreciated. This is 
why here we concentrate on pulse-coupled oscillators. 

Finally, for what concerns the topology of the interactions, 
it is known that they can heavily influence the dynamics of the 
neural systems leading to the emergence of new collective phe- 
nomena even in weakly connected networks (Timme, 2006), or 
of various types of chaotic behavior, ranging from weak chaos 
for diluted systems (Popovych et al., 2005; Olmi et al., 2010) 
to extensive chaos in sparsely connected ones (Monteforte and 
Wolf, 2010; Luccioli et al., 2012). We will, however, limit our 
analysis to globally coupled identical oscillators, which provide 
a much simplified, but already challenging, test bed. The high 
symmetry of the corresponding evolution equations simplifies the 
identification of the stationary solutions and the analysis of their 
stability properties. The two most symmetric solutions are: (1) 
the fully synchronous state, where all oscillators follow exactly the 
same trajectory; (2) the splay state (also known as "ponies on a 
merry-go-round," antiphase state or rotating waves) (Hadley and 
Beasley, 1987; Ashwin et al, 1990; Aronson et al, 1991), where 
the oscillators still follow the same periodic trajectory, but with 
different (evenly distributed) time shifts. The former solution is 
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the simplest representative of the broad class of clustered states 
(Golomb and Rinzel, 1994), where several oscillators behave in 
the same way, while the latter is the prototype of asynchronous 
states, characterized by a smooth distribution of phases (Renart 
etal, 2010). 

In spite of the many restrictions on the mathematical setup, 
the stability of the synchronous and splay states still depend 
significantly on additional features such as the synaptic response- 
function, the velocity field, and the presence of delay in the pulse 
transmission. As a result, one can encounter splay states that 
are either strongly stable along all directions, or that present 
many almost-marginal directions, or, finally, that are marginally 
stable along various directions (Nichols and Wiesenfield, 1992; 
Watanabe and Strogatz, 1994). Several analytic results have been 
obtained in specific cases, but a global picture is still missing: the 
goal of this paper is to recompose the puzzle, by exploring the 
role of the velocity field (or, equivalently, of the phase response 
curve) and of the shape of the transmitted post-synaptic poten- 
tials. Although we are neither going to discuss the role of delay 
nor that of the network topology, it is useful to recall the stabil- 
ity analysis of the synchronous state in the presence of delayed 
8-pulses and for arbitrary topology, performed by Timme and 
Wolf in Timme and Wolf (2008). There, the authors show that 
even the complete knowledge of the spectrum of the linear oper- 
ator does not suffice to address the stability of the synchronized 
state. 

The stability analysis of the fully synchronous regime is far 
from being trivial even for a globally coupled network of oscil- 
lators with no delay in the pulse transmission: in fact, the pulse 
emission introduces a discontinuity which requires separating the 
evolution before and after such event. Moreover, when many neu- 
rons spike at the same time, the length of some interspike intervals 
is virtually zero but cannot be neglected in the mathematical 
analysis. In fact, the first study of this problem was restricted to 
excitatory coupling and 8-pulses (Mirollo and Strogatz, 1990). In 
that context, the stability of the synchronous state follows from 
the fact that when the phases of two oscillators are sufficiently 
close to one another, they are instantaneously reset to the same 
value (as a result of a non-physical lack of invertibility of the 
dynamics). The first, truly linear stability analyses have been per- 
formed later, first in the case of two oscillators (van Vreeswijk 
et al., 1994; Hansel et al., 1995) and then considering 8-pulses 
with continuous PRCs (Goel and Ermentrout, 2002). Here, we 
extend the analysis to generic pulse-shapes and discontinuous 
PRCs [such as for leaky integrate and fire (LIF) neurons]. 

As for the splay states, their stability can be assessed in two 
ways: (1) by assuming that the number of oscillators is infi- 
nite (i.e., taking the so called thermodynamic limit) and thereby 
studying the evolution of the distribution of the membrane 
potentials — this approach is somehow equivalent to dealing with 
(macroscopic) Liouville-type equations in statistical mechanics; 
(2) by dealing with the (microscopic) equations of motion for 
a large but finite number N of oscillators. As shown in some 
pioneering works (Kuramoto, 1991; Treves, 1993), the former 
approach corresponds to develop a mean field theory. The result- 
ing equations have been first solved in Abbott and van Vreeswijk 
(1993) for pulses composed of two exponential functions, in the 



limit of a small effective coupling [A small effective coupling can 
arise also when PRC has a very weak dependence on the phase 
(see section 3)]. Here, following Abbott and van Vreeswijk (1993), 
we extend the analysis to generic pulse-shapes, finding that sub- 
stantial differences exist among 8, exponential and the so-called 
a-pulses (see the next section for a proper definition). 

Direct numerical studies of the linear stability of finite net- 
works suggest that the eigenfunctions of the (Floquet) operator 
can be classified according to their wavelength I (where I refers 
to the neuronal phase — see section 4.1 for a precise definition). In 
finite systems, it is convenient to distinguish between long (LW) 
and short (SW) wavelengths. Upon considering that I = n/N 
(1 < n < N), LW can be identified as those for which n N, 
while SW correspond to larger n values. Numerical simulations 
suggest also that the time scale of a LW perturbation typically 
increases upon increasing its wavelength, starting from a few mil- 
liseconds (for small n values) up to much longer values (when n 
is on the order of the network size N) which depend on "details" 
such as the continuity of the velocity field, or the pulse shape. On 
the other hand, SW are characterized by a slow size-dependent 
dynamics. 

For instance, in LIF neurons coupled via a-pulses, it has been 
found (Calamai et al., 2009) that the Floquet exponents of LW 
decrease as 1 /I 2 (for large I), while the time scale of the SW 
component is on the order of N 2 . In practice the LW spectral 
component as determined from the finite N analysis coincides 
with that one obtained with the mean field approach (i.e., tak- 
ing first the thermodynamic limit). As for the SW component, it 
cannot be quantitatively determined by the mean-field approach, 
but it is nevertheless possible to infer the correct order of mag- 
nitude of this time scale. In fact, upon combining the 1 /I 2 decay 
(predicted by the mean-field approach) with the observation that 
the minimal wavelength is 1 /N, it naturally follows that the SW 
time scale is N 2 , as analytically proved in Olmi et al. (2012). 
Furthermore, it has been found that the two spectral components 
smoothly connect to each other and the predictions of the two 
theoretical approaches coincide in the crossover region. 

It is therefore important to investigate whether the same agree- 
ment extends to more generic pulse shapes and velocity fields. The 
finite-N approach can, in principle, be generalized to arbitrary 
shapes, but the analytic calculations would be quite lengthy, due 
to the need of distinguishing between fast and slow scales and the 
need of accounting for higher order terms. For this reason, here 
we limit ourselves to give a positive answer to this question with 
the help of numerical studies. 

The only, important, exception to this scenario is obtained for 
quasi 8-like pulses (Zillmer et al, 2007), i.e., for pulses whose 
width is smaller than the average time separation between any two 
consecutive spikes, in which case all the SW eigenvalues remain 
finite for increasing N. 

In section 2 we introduce the model and derive the corre- 
sponding event-driven map, a necessary step before undertaking 
the analytic calculations. Section 3 is devoted to a perturbative 
stability analysis of the splay state in the infinite-size limit for 
generic velocity fields and pulse shapes. The following section 4 
reports a discussion of the stability in finite networks. There we 
briefly recall the main results obtained in Olmi et al. (2012) for 
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the splay state and we extensively discuss the method to quantify 
the stability of the fully synchronous regime. The following two 
sections are devoted to a numerical analysis of various setups. 
In section 5 we study splay states in finite networks for generic 
velocity fields and three different classes of of pulses, namely, 
with finite, vanishing (sal/JV), and zero width. In section 6 we 
study periodically forced networks. Such studies show that the 
scaling relations derived for the splay states apply also to such a 
microscopically quasi-periodic regime. A brief summary of the 
main results together with a recapitulation of the open problem 
is finally presented in section 7. In the first appendix we derive 
the Fourier components needed to assess the stability of a splay 
state for a generic PRC. In the second appendix the evaporation 
exponent is determined for the synchronous state in LIF neurons. 

2. THE MODEL 

The general setup considered in this paper is a network of N 
identical pulse-coupled neurons (rotators), whose evolution is 
described by the equation 

X? =FQd)+gE(t), j=l,...,N (1) 

where X } represents the membrane potential, g is the coupling 
constant and the mean field E(t) denotes to the synaptic input, 
common to all neurons in the network. When X 1 reaches the 
threshold value X? = 1, it is reset to X? = 0 and a spike contributes 
to the mean field £ in a way that is described here below. The 
resetting procedure is an approximation of the discharge mecha- 
nism operating in real neurons. The function F(X) (the velocity 
field) is assumed to be everywhere positive, thus ensuring that 
the neuron is repetitively firing. For Fq(X) = a — X the model 
reduces to the well-known case of LIF neurons. 

The mean field E(t) arises from the linear superposition of the 
pulses emitted by the single neurons. In order to describe its time 
evolution, it is sufficient to introduce a suitable ordinary differen- 
tial equation (ODE), such that its Green function reproduces the 
expected pulse shape, 

B ( «=|J + | £ Ht-t n ), (2) 

i n\t n <t 

where the superscript (i) denotes the ith time derivative, L the 
order of the differential equation and K = Y[, ( — <*> being the 
poles of the differential equation), so as to ensure that the sin- 
gle pulses have unit area (for N = 1). The 8-functions appearing 
on the right hand side of Equation (2) correspond to the spikes 
emitted at times {f„}: each time a spike is emitted, the term Ef- L ~^ 
has a finite jump of amplitude K/N. Therefore L controls the 
smoothness of the pulses: L — 1 is the order of the lowest deriva- 
tive that is discontinuous. L = 0 corresponds to the extreme case 
of 8-pulses with no field dynamics; L = 1 corresponds to discon- 
tinuous exponential pulses; L = 2 (with ai = 0.2) to the so-called 
a-pulses (E s (t) = a 2 te~ at ). Since a-pulses will be often referred 
to, it is worth being a little more specific. In this case, Equation (2) 
reduces to 

a 2 

E(t) + 2aE(t) + a 2 E(t) = — 8 ( f - f «)> ( 3 ) 

n\t„ < t 



and it is convenient to transform this equation into a system of 
two ODEs, namely 

a 2 

E = P-aE, P + aP=— s ( f - f «). ( 4 ) 

n\t n <t 

where we have introduced, for the sake of simplicity, the auxiliary 
variable P = aE + E. 

2.1. EVENT-DRIVEN MAP 

By following Zillmer et al. (2006) and Calamai et al. (2009), it is 
convenient to pass from a continuous — to a discrete-time evolu- 
tion rule, by deriving the event-driven map which connects the 
network configuration at consecutive spike times. For the sake 
of simplicity, in the following part of this section we refer to 
a-pulses, but there is no conceptual limitation in extending the 
approach to L > 2. 

By integrating Equation (4), we obtain 

E n+1 = E n e- aT " + P n T n e- aT " (5) 
a 2 

P n + 1 =P n e- aT " + -, (6) 

where we have taken into account the effect of the incoming 
pulse (see the term a 2 /N in the second equation) while T n = 
t n + 1 — t n is the interspike interval; t n+ 1 corresponds to the time 
when the neuron with the largest membrane potential reaches the 
threshold. 

Since all neurons follow the same first-order differential equa- 
tion (this is a mean-field model), the ordering of their membrane 
potentials is preserved [neurons "rotate" around the circle [0, 1] 
without overtaking each other (Jin, 2002)]. It is, therefore, con- 
venient to order the potentials from the largest to the smallest 
one and to introduce a co-moving reference frame, i.e., to shift 
backward the label;', each time a neuron reaches the threshold. By 
formally integrating Equation (1), 

1 i+ 1 ^ e~ r » - e- aT » ( P n \ 

+ a— 1 \ a — 1 / 



Moreover, since X\ is always the largest potential, the interspike 
interval is defined by the threshold condition 

X l n {T n ,E n ,P n ) = l. (8) 

Altogether, the model now reads as a discrete-time map, involv- 
ing N + 1 variables, E„, P n , andX^, (1 < j < N), since one degree 
of freedom has been eliminated as a result of having taken the 
Poincare section (X^ = 0 due to the resetting mechanism). The 
advantage of the map description is that we do not have to deal 
any longer with 8 -like discontinuities, or with formally infinite 
sequences of past events. 

In this framework, the splay state is a fixed point of the event- 
driven map. Its coordinates can be determined in the following 
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way. From Equation (5), one can express P and £ as a function of 
the yet unknown interspike interval T, 



N 



(1 



e- aT r l 



: TP(e' 



1)" 



(9) 



The value of the membrane potentials X are then obtained by 
iterating backward in; Equation (7) (the n dependence is dropped 
for the fixed point) starting from the initial condition X N = 0. 
The interspike interval T is finally obtained by imposing the con- 
dition X° = 0. In practice the computational difficulty amounts 
to finding the zero of a one dimensional function and, even 
though T{X) + 1 , T) can, in most cases, be obtained only through 
numerical integration, the final error can be very well kept under 
control. 

3. THEORY (N = oo) 

The stability of a dynamical state can be assessed by either first 
taking the infinite-time limit and then the thermodynamic limit, 
or vice versa. In general it is not obvious whether the two methods 
yield the same result and this is particularly crucial for the splay 
state, as many eigenvalues tend to 0 for N — >■ oo. In this section 
we discuss the scenarios that have to be expected when the ther- 
modynamic limit is taken first. We do that by following Abbott 
and van Vreeswijk (1993). 

As a first step, it is convenient to introduce the phase-like 
variable 



» = f 

Jo 



dx 

GOO' 



0 < y' < 1 



(10) 



where, for later convenience, we have defined G(X) = g + 
TqF(X), To = NT being the period of the splay state (i.e., the 
single-neuron interspike interval). The phase y' evolves according 
to the equation 



dy' 
dt 



E + 



gs(0 
G(X(f)) 



(11) 



where E = 1 /To is the amplitude of the field in the splay state, 
e(t) = E(t) — E. In the splay state, since e = 0, y' grows lin- 
early in time, as indeed expected for a well-defined phase. In the 
thermodynamic limit, the evolution is ruled by the continuity 
equation 



dt 



97 

dy 



where p(y, t)dy is the fraction of neurons whose phase y' lies in 
(y, y + dy) at time f, and 



J(y, f) 



E + 



G(X(y))J 



p(y, t) 



(13) 



is the corresponding flux. As the resetting implies that the out- 
going flux /(l, f) (which coincides with the firing rate) equals 



complemented with the boundary condition /(0, f) = J(l, t). 
Finally, in this macroscopic representation, the field equation 
writes 



L-l 



= J2 fl«e (i) + K{J{1, t)-E), 



(14) 



while the splay state corresponds to the fixed point p = 1, g = 0, 
J = E. The smoothness of the splay state justifies the use of a par- 
tial differential equation such as (Equation 12). Its stability can be 
studied by introducing the perturbation j(y, t) 

j(y,t)=J{y,t)-E, (15) 

and linearizing the continuity equation, 

-± = — 5 16 

dt G(X(y)) dt dy 

while the field equation simplifies to 

i-i 

e (i) =^ai£® +Kj(l,f). (17) 

By now introducing the Ansatz 

j(y, t) = ;/(y)e If e(y, t) = e/(y)e Xr , (18) 

in Equations (16) and (17) and, thereby solving the resulting 
ODE, one can obtain an implicit expression for jfiy), 



jf{y) = e-^ S 



1 + 



gKX j f (l) 



Ei\ k= i^ + a k)Jo G(X(z)) 



dz 



\z/E 



where — ak and K are defined as below Equation (2). By imposing 
the boundary condition for the flux, jf(l) = jf(0) = 1, one finally 
obtains the eigenvalue equation (Abbott and van Vreeswijk, 
1993), 

(e^-.)ri(H,) = ef^. (19) 
V J Li E ^ G(X(y)) 



(12) I n me case °f a constant G(X(y)) = a, L eigenvalues correspond 
to the zeroes of the following polynomial equation 



n ^ + ^ 



(20) 



For g = 0 such solutions are the poles which define the field 
dynamics, while for g = a, X = 0 is a solution: this corresponds 
to the maximal value of the (positive) coupling strength beyond 
which the model does no longer support stationary states, as 



the incoming flux at the origin, the above equation has to be the feedback induces an unbounded growth of the spiking rate. 
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Besides such L solution, the spectrum is composed of an infinite 
set of purely imaginary eigenvalues, 



2:tin 

"k = 2itin£ = n^O. 



(21) 



The existence of such marginally stable directions reflects the fact 
that all y' phases experience the same velocity field, indepen- 
dently of their current value ( see Equation 1 1 ) , so that no effective 
interaction is present among the oscillators. In the limit of small 
variations of G(X(y)), one can develop a perturbative approach. 
Here below, we proceed under the more restrictive assumption 
that the coupling constant g is itself small: we have checked that 
this restriction does not change the substance of our conclusions, 
while requiring a simpler algebra. 

A small g value implies that "k is close to 2jtmE and thereby 
expand the exponential in Equation (19). Up to first order, we 
find 



X„ = 2jtmE 



1 + 



gK(A n + iB n ) 
l\l =l (2nmE + a k ) 



where 



(A n + iB„) 



I 



Jlstny 



dy 

o G(X(y)) 



(22) 



(23) 



are the Fourier components of the phase-response curve 
l/G(X(y)). 

In order to estimate the leading terms of the real part of X„ in 
the large n limit, let us rewrite Equation (22) as 



K = iy n + gKy„ B " + lA " FT ( ak _ (24) 

where y n = litnE = (27tn)/To. Since y„ is proportional to n, the 
leading terms in the product at numerator of Equation (24) are 



k=\ 



1-1 



(25) 



where S = Ylk= i ot* while the leading term in the product at 
denominator in Equation (24) is }/ 2i . Accordingly, the main con- 
tribution to the real part of the eigenvalues is, in the case of 
even L, 



Re[k n ] ~ gK(-\) 1 ' 2 
and, for odd L, 

Re{\ n ) ~^(-l)( i+3 >/ 2 



SA n B n 

' n Yn 



A n SB n 
i n . 



Yn 



(26) 



(27) 



An exact expression for the Fourier components A n and B n 
appearing in Equation (23) can be derived in the large n limit. 



In particular, the integral over the interval [0, 1] appearing in 
Equation (23) can be rewritten as a sum of integrals, each 
performed on a sub-interval of vanishingly small length l/n. 
Furthermore, since the phase-response 1/G has a limited varia- 
tion within each sub-interval, it can be replaced by its polynomial 
expansion up to second order. Finally, as shown in Appendix A, 
the following expression are obtained at the leading order in 1 /n 
for a discontinuous F(X) 



An 



-T, 



II 



4it 2 n 2 

To 
2im 



F(l) f(0) 
G(l) 2 ~ G(0) 2 

f(i)-f(Q) 

G(1)G(0) 



Therefore, for even L, the leading term for n — > oo is 



Re{-k„} 



1)1/2 (F( 0)_ F (1)) 



(2jt«) L G(l)G(0) 



(28) 
(29) 

(30) 



For even L, the stability of the short-wavelength modes (large n) 
is controlled by the sign of (-F(O) — F(l)): for even (odd) 1/2 and 
excitatory coupling, i.e., g > 0, the splay state is stable whenever 
-F(l) > F(0) (-F(l) < f(0)). Obviously the stability is reversed for 
inhibitory coupling. 

Notice that for L = 0, i.e., 8-spikes, the eigenvalues do not 
decrease with n, as previously observed in Zillmer et al. (2007). 
This is the only case where all modes exhibit a finite stability even 
in the thermodynamic limit. 

For odd L, the real part of the eigenvalues is 



Re{\„} 



*KT£(-1) 



(L+D/2 



(2TtM)( i+1 > 

\ F(l) F'(0) 



(31) 



1 G(l) 2 



G(0)2 



sr 0 



F{1) - f (0) 
G(1)G(0) 



in this case the value of F(X) and of its derivative F'(X) at the 
extrema mix up in a non-trivial way. 

Finally, as for the scaling behavior of the leading terms we 
observe that 



Re{\„} 



L+l 

2 



(32) 



where |_-J stays for the integer part of the number. Therefore the 
scaling of the short-wavelength modes for discontinuous F(X) is 
dictated by the post-synaptic pulse profile. 

For a continuous but non-differentiable F(X), (i.e., -F'(l) ^ 
-F'(O)), if L is even, it is necessary to go two orders beyond in the 
estimate of the Fourier coefficients (see Appendix A). As a result, 
the eigenvalues scale as 



Re{X n } a n 



-(1+2) 



(33) 



For odd L, it is instead sufficient to assume F(0) = F(l) in 
Equation (31). 

Altogether, we have seen that the non-smoothness of both 
the post-synaptic pulse and of the velocity field (or, equivalently, 
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of the phase response curve) play a crucial role in determining 
the degree of stability of the splay state. The smoother are such 
functions and the slower short-wavelength perturbations decay, 
although the changes occur in steps which depend on the parity 
of the order of the discontinuity (at least for the pulse struc- 
ture). Moreover, the overall stability of the spectral components 
depends in a complicate way on the sign of the discontinuity itself. 

4. THEORY (FINITE N) 
4.1. THE SPLAY STATE 

The stability for finite N can be investigated by linearizing 
Equations (5-7). A thorough analysis has been developed in Olmi 
et al. (2012); here we limit ourselves to review the key ideas as a 
guide for the numerical analysis. 

We start by introducing the vector W = ({x 1 }, e, p) 
(j=l,N—l), whose components represent the infinites- 
imal perturbations of the solution {X 1 }, E, P. The Floquet 
spectrum can be determined by constructing the matrix A which 
maps the initial vector W(0) onto W(T), 



W(T) = AW(0) 



(34) 



been found that the results reported in Abbott and van Vreeswijk 
(1993) match does obtained for 1 < k < N in Olmi et al. (2012). 
The latter one corresponds to the SW component: it depends on 
the system size and cannot, indeed, be derived from the mean field 
approach discussed in the previous section. In the next section, 
we illustrate some examples that go beyond the analytic studies 
carried out in Olmi et al. (2012). 

4.2. THE SYNCHRONIZED STATE 

In this section we address the problem of measuring the stability 
of the fully synchronized state for a generic oscillator dynamics 
F(x). The task is non-trivial, because of the resetting mecha- 
nism, which acts simultaneously on all neurons. On the one side, 
we extend the results obtained in Goel and Ermentrout (2002) 
which are restricted to a continuous PRC, on the other side we 
extend the results of Mirollo and Strogatz (1990) which refer to 
excitatory coupling and 8 pulses. In order to make the analysis 
easier to understand we start considering a-pulses. Other cases 
are discussed afterward. 

The starting point amounts to writing the event driven map in 
a comoving frame, 



where T corresponds to the time separation between two con- 


K+\ 




A„ 


1 


secutive spikes. This is done in two steps, the first of which 




-aT„ 




corresponds to evolving the components of a Cartesian basis 


E n + l 


= E n c 


+ 


according to the equations obtained from the linearization of 










Equations (1, 4) (in the comoving reference frame), 


Pn+l 


= P n e 




+ 



V 



dF 



axj+i 
€ = p — ote, p 



x j+1 +ge, j = 2,...,N ^ = 0 



-ap. 



(35) 



The second step consists in accounting for the spike emission, 
which amounts to add the vector 



U=[{X>(T)},E(T),P(T)]x 



(36) 



where x is obtained from the linearization of the threshold 
condition (8), 



1 



3X 1 dX 1 

€ -I p [ . 

dE dP ) X 1 



(37) 



The diagonalization of the resulting matrix A, gives N + 1 
Floquet eigenvalues nt, which we express as 



e '<l>k e To(y*k + i«>k)/N 



(38) 



where <pt = ^jf, k = 1, . . . , N — 1, and </>jv = </>n-i = 0, while 
Xk and tojt are the real and imaginary parts of the Floquet expo- 
nents. The variable 0^ plays the role of the wavenumber k in the 
linear stability analysis of spatially extended systems. 

Previous studies (Olmi et al., 2012) have shown that the spec- 
trum can be decomposed into two components: (1) k ~ 0(1); 
(2) k/N ~ 0(1). The former one is the LW component and can 
be directly obtained in the thermodynamic limit (see the previ- 
ous section). For L = 2 and ai = ct2 (i-e-> for a pulses), it has 



N 



(39) 
(40) 
(41) 



where the function T is obtained by formally integrating the 
equations of motion over the time interval T n . Notice that the 
field dynamics has been, instead, explicitly obtained from the 
exact integration of the equations of motion [compare with 
Equations (3, 4)]. The interspike time interval T n is finally deter- 
mined by solving the implicit equation 



F{X l n ,E n ,P n ,T n ) = \. 



(42) 



In order to determine the stability of the synchronized state, it is 
necessary to assume that the neurons have an infinitesimally dif- 
ferent membrane potentials, even though they coincide with one 
another. As a result, the full period must be broken into N steps. 
In the first one, of length T, all neurons start in X = 0 and arrive 
at 1, but only the "first" reaches the threshold; in the following 
N — 1 steps, of 0-length, one neuron after the other passes the 
threshold and it is accordingly reset in 0. 

With this scheme in mind we proceed to linearize the equa- 
tions, writing the evolution equations for the infinitesimal per- 
turbations x'n, 6„, p n , and x n around the synchronous solution. 
From Equations (39-41) we obtain, 

4+ 1 = Fx(j + D*!, + ' + FeQ + + 

TA] + l)pn + Ft(] + l)t« 1 < j < N (43) 

6„+i = e- aT t „ + Te- aT p n - 

(a~E - P n e- aT ) X„ (44) 

Pn+ i = e- aT p„ - aP n e- aT x n . (45) 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



February 2014 | Volume 8 | Article 8 | 6 



Olmi et al. 



Linear stability in networks of pulse-coupled neurons 



with the boundary condition 1 = 0 (due to the reset mech- 
anism) and where the subscripts X, E, P, and T denote a partial 
derivative with respect to the given variable. Moreover, the depen- 
dence on j + 1 is a shorthand notation to remind that the various 
derivatives depend on the membrane potential of the + l)st 
neuron. Finally, we have left the n-dependence in the variable 
P as it changes (in a 2 /N steps, when the neurons progressively 
cross the threshold), while E refers to the field amplitude, which, 
instead, stays constant. 

The above equations must be complemented by the condition 



-T x x l n +T E e n +Tpp n , 



(46) 



where T z = FzW/FtO) (Z = X, E, P). Equation (46) is 
obtained by differentiating Equation (42) which defines the 
period of the splay state. 

We now proceed to build the Jacobian for each of the N steps, 
starting from the first one. In order not to overload the notations, 
from now on, the time index n corresponds to the step of the pro- 
cedure. It is convenient to order all the variables, starting from x 1 
( j = 1, N — 1), and then including 6 andp, into a single vector, so 
that the evolution is described by an (N + 1) x (N + 1) matrix 
with the following structure, 



7V(n) 



T(n) 0 
*(«) fi(n) 



(47) 



where 0 is an (N — 1) x 2 null matrix; T(n) is a quadratic 
(N — 1) x (N — 1) matrix, whose only non-zero elements are 
those in the first column and along the supradiagonal; is 
a 2 x (N — 1) matrix whose elements are all zero except for the 
first column; finally £2 (n) is a 2 x 2 matrix. 

Since in the first step all neurons start from the same position 
X = 0, one can drop the j dependence in T . With the help of 
Equations (46, 43) 



r(i) ;il = -t x 
r(% +1 = Tx 

Moreover, with the help of Equations (44-46) 



(48) 



*(l)n = -(aE-Pe- aT ^T x 
*(1) 12 = -aPe- aT T x , 



(49) 



where we have also made use that Pi = P. Finally, 
Q(\) u = e - aT - (aE - Pe- aT ^j T E , 

Q(l) u = Te- aT - (aE - P e - aT ^j % 
h„-aTr 



a(\) 21 =-aPe- al T E , 
J2(l) 22 = e - aT - aPe- aT T P . 



(50) 



In the next steps, T n vanishes, so that Te = J~p = 0, while Tx = 
1 and JFt(1) = ^CO + gE '■= V 1 . Moreover, Ttif) depends on 



whether the jth neuron has passed the threshold or not. In the for- 
mer case .FtO + 1) = F(0) + gE := Vb> otherwise ZFrij + 1) = 
V 1 . As a result, 



r (»);,! = -v'/v 1 

r («)/,;+ 1 = i 



(51) 



where V = V° if j < n and V = V 1 , otherwise. At the same 
time, from the equations for the field variables, we find that 



*(«)n 



*(«)i 



a ~E-(P+(n-l)^) 
V 1 

ct(P+(H-l)g) 
V 1 



(52) 



while Q. (n) reduces to the identity matrix. 

From the multiplication of all matrices, we find that the 
structure is preserved, namely 



Af(N) ■ ■ ■ N(2)Af(l) 



A 0 

* £2(1) 



(53) 



where *(h) is a 2 x (N — 1) matrix, whose elements are all zero 
except for those of the first column, namely 

*n = *(l)n + *(«)n 

*12 = *(l)l2 + *(«)l2 

Furthermore, A is a diagonal matrix, with 



A„ 



V° F{0)+gE 
■ SX7~i- = '. ~ ex P 



V 1 



F(l)+gE 



f 

Jo 



dtF'(X(t)) 



(54) 



Therefore, it is evident that the stability of the orbit is measured 
by the diagonal elements A« together with the eigenvalues of Q. 
which are associated to the pulse structure. In practice, Tx cor- 
responds to the expansion rate from X = 0 to X = 1 under the 
action of the mean field E and we recover a standard result in 
globally coupled identical oscillators: the spectrum is degenerate, 
all eigenvalues being equal and independent of the network size. 
The result is, however, not obvious in this context, due to the care 
that is needed in taking into account the various discontinuities. 
We have separately verified that the same conclusion holds for 
exponential spikes. 

The stability of the synchronized state can be also addressed by 
determining the evaporation exponent A e (van Vreeswijk, 1996; 
Pikovsky et al, 2001), which measures the stability of a probe 
neuron subject to the mean field generated by the synchronous 
neurons with no feedback toward them. By implementing this 
approach for a negative perturbation, van Vreeswijk found that 
A e is equal to A« (for a-functions). By further assuming that 
F' < 0, he was able to prove that the synchronized state is sta- 
ble for inhibitory coupling and sufficiently small a-values. The 
situation is more delicate for exponential pulse-shapes. As shown 
in di Volo et al. (2013), A e > 0 (A e < 0) depending whether the 
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perturbation is positive (negative). In this case, the Floquet expo- 
nent reported in Equation (54) coincides with the evaporation 
exponent estimated for negative perturbations. In Appendix B. 
we show that the difference between the left and right stability is to 
be attributed to the discontinuous shape of the pulse: no anomaly 
is expected for a pulses. 

5. NUMERICAL ANALYSIS 

The theoretical approaches discussed in the previous sections 
allow determining: (1) the SW components of the Floquet spec- 
trum for discontinuous velocity fields; (2) the leading LW expo- 
nents directly in the thermodynamic limit for generic velocity 
fields and pulse shapes, in the weak coupling limit. It would be 
possible to extend the finite N results to other setups, but we do 
not think that the effort is worth, given the huge amount of tech- 
nicalities. We thus prefer to illustrate the expected behavior with 
the help of some simulations which, incidentally, cover a wider 
range than possibly accessible to the analytics. 

More precisely, in this and the following section we study the 
models listed in Table 1 in a standard set up (splay states) and 
under the effect of periodic external perturbations. 

5.1. FINITE PULSE WIDTH 

Here, we discuss the stability of the splay state for different degrees 
of smoothness of the velocity field at the borders of the unit 
interval for post-synaptic pulses of a-function type. 

We start from discontinuous velocity fields. They have been the 
subject of an analytic study which proved that the SW component 
scales as 1/AT 2 (Olmi et al., 2012). The data reported in Figure 1A 
for F\ (X) confirms the expected scaling: the agreement with the 
theoretical curve derived in Olmi et al. (2012) is impressive over 
the entire spectral range, while the mean field Equation (30) gives 
a very good estimation of the spectrum except for the shortest 
wavelengths, where it overestimates the numerical data. The mean 
field approximation turns out to be more accurate for continuous 
velocity fields (with a discontinuity of the first derivative at the 



Table 1 | In the first column is reported the list of the velocity fields 
F\X) analyzed in the paper. All the considered fields are everywhere 
positive within the definition interval Xe[0,1], thus ensuring that the 
neuron is supra-threshold. The second column refers to the 
continuity properties of the fields within the interval [0,1]. 



Velocity field Continuity properties 



F 0 (X) = 


a — 


X 


Discontinuous 


Fi(X) = 


a - 


X(X-0.7) 


Discontinuous 


F 2 (X) = 


a - 


0.25sinCrtX) 


c<°> 


F 3 (X) = 


a + 


X(X- 1) 


C (0) 


F 4 (X) = 


a - 


0.25sin(uX)cos 2 0itX) 


c<°> 


FsVO = 


a - 


0.25sin(2nX) cos 2 (2^X) 


C (oo) 


F 6 (X) = 


a - 


0.25sin(2nX)e cos(2,tX) 


C (oo) 


FAX) = 


a — 


1 _|_ e 2sin(27iX) 


C (oo) 



The function is labeled as discontinuous if F(0) ^ F(1); it is C l0> if F(0l = F(l) but 
F'(0)£ F'iV andC lv ifF(0)= F(V andF'(0)= F'(1). F(X)isC ,ccl if it is infinitely 
differentiable and each derivative is continuous at the extrema of the definition 
interval. 



borders of the definition interval). Indeed the agreement between 
the theoretical expression Equation (A10) and the numerical data 
is very good for the entire range [see Figure IB which refers to 
MX)]. 

The numerical Floquet spectra for fields that are C'°\ but 
not C (1) (F(Q) = F(l), F'(0) / -F'(l)), are reported in Figure 2 
[the curves in panels (A, B) refer to Fi(X) an d F4, respectively]. 
For these velocity fields, we have also verified that the spectra 
scale as 1/N 4 , confirming the observation reported in Calamai 
et al. (2009) for a different velocity field with the same analyti- 
cal properties. The data displayed in Figures 2A,B refer to the LW 
components: they indeed confirm to be independent of the sys- 
tem size and scale as 1 /k 4 (see the dashed line) as predicted by the 
perturbative theory discussed in section 3. 

The spectra reported in the other two panels refer to analytic 
velocity fields: in all cases the initial part of the Floquet spectra is 
again independent of N and scales approximately exponentially 
with k, confirming that the scaling behavior of the exponents 
is related to the analyticity of the velocity field. The fluctuating 
background with approximate height 10~ 12 is just a consequence 
of the finite numerical accuracy. This is the reason why we did not 
dare to estimate the SW components that would be exceedingly 
small. 

5.2. VANISHING PULSE-WIDTH 

Here, we analyze the intermediate case between finite pulse-width 
and S-like impulses. Similarly to what done in Zillmer et al. 
(2007) for the LIF, we consider a pulses, where a = fW, with f5 
independent of N. 

In Figure 3A we report the spectra for a discontinuous veloc- 
ity field, F\(x). In this case the Floquet spectra remain finite, so 
that the corresponding states remain robustly stable even in the 
thermodynamic limit. Also in this case the agreement with the 
theoretical expression reported in Equation (7) in Olmi et al. 
(2012) is extremely good, while Equation (30) overestimates the 
spectra for large phases. The field considered in panel (b) (-F2PO) 
is but notC' 1 '. In this case, the Floquet spectra scale as l/N: 
this scaling is predicted by the analysis reported in section 3 and 
the whole spectrum is very well reproduced by Equation (A10). 




FIGURE 1 | Floquet spectra for a-pulses for (A) a discontinuous field 
Fi(X) and (B) a continuous field Ft(X). The orange dotted line in (A) 
represents the theoretical curve estimated by using Equation (7) in Olmi 
et al. (2012), while the dashed maroon curve represents the theoretical 
curve estimated by using Equation (30) in section 3. In (B) the dashed 
maroon curve is calculated by using Equation (A10). All data refer to a= 1.3 
and a = 3. 
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1000 



FIGURE 2 | Floquet spectra for a-pulses for two continuous sinusoidal 
fields, namely F 2 [X) (A) and F 4 (X) (B); and two analytic fields, namely 
F${X) (C) and Fj[X) (D). The dashed blue line in (B) indicates a scaling 
1 /k 4 . All data refer to a = 1 .3 and a = 3. 




FIGURE 3 | Floquet spectra for (i-pulses with a discontinuous field 
[F,(X)] (A) and a C m field [F 2 (X)] (B). The orange dotted line in (A) 
represents the theoretical curve estimated by using Equation (7) in Olmi 
et al. (2012). The dashed line in (A) [resp. (B)] represents the theoretical 
curve computed by using Equation (30) [resp. Equation (A10)] for fi-pulses. 
The data refer to a = 1 .3 and p = 0.03. 



Last but not least, we have studied an analytic field, namely 
Fj{X). In this case the Floquet spectra appear to scale exponen- 
tially to zero with the wavevector k, similarly to what observed for 
the finite pulse width, as shown in Figure 4. 

5.3. f> PULSES 

Finally we considered the case of 8-pulses: whenever the potential 
X reaches the threshold value, it is reset to zero and a spike is sent 
to and instantaneously received by all neurons. We studied just 
two cases: (1) the analytic field Fj(X); (2) a leaky integrate-and 
fire neuron model with Fq(X). The results, obtained for inhibitory 
coupling [since the splay state is known to be stable only in such 
a case (van Vreeswijk, 1996; Zillmer et al, 2006)] are consistent 
with the expectation for the f5 model. 

In particular we found, in the analytic case ( 1 ), that the Floquet 
spectra decay exponentially to zero. The exponential scaling is not 
altered if a phase shift i; is introduced in the velocity field (i.e., for 
F(X) = a- 1 + e 2sin(2nX+ V). In the case of the LIF model (F 0 ), 




hid 



FIGURE 4 | Floquet spectra for ({-pulses for the analytic field Fj[X). The 

data refer to a = 1 .3 and p = 0.03. 



we already know that the Lyapunov spectrum tends, in the 8-pulse 
limit, to Zillmer et al. (2007) 



lim X,t 



1 , / a 
-l+-ln — 
To \a- 



(55) 



This result is confirmed by our simulations which also reveal that 
the splay state is stable even for small, excitatory coupling values, 
extending previous results limited to inhibitory coupling (Zillmer 
etal, 2006). 

6. PERIODIC FORCING 

In this section we numerically investigate the scaling behavior of 
the Floquet spectrum in the presence of a periodic forcing, to test 
the validity of the previous analysis in a more general context. We 
have restricted our studies to splay-state-like regimes, where it is 
important to predict the behavior of the many almost marginally 
stable directions. Moreover, we have considered only the smooth 
a-pulses. In this case, the dynamical equations read 

X = F(X) + gE + A cos((p), j=l,...,N, 
E = P-aE, (56) 
P = -aP, 

(p = (JO. 

They have been written in an autonomous form, since it is more 
convenient to perform the Poincare section according to the 
spiking times, rather than introducing a stroboscopic map. The 
interspike interval is determined by the equation 



T: 



/' 



dX l 



F(X l )+gE + Acos(q>)' 



(57) 



where X 1 is the membrane potential of the first neuron (the 
closest to threshold), and X 0 ia is its initial value. 

We analyzed only those setups where the unperturbed splay 
state is stable. More precisely: the two discontinuous fields Prj (X) 
and Fi (X) , the two C (0) fields (F 2 (X) and F 3 (X) ) , and the analytic 
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FIGURE 5 | Lyapunov spectra for neurons forced by an external 
periodic signal, we observe the scaling 1/A/ 2 for the discontinuous 
velocity fields (A) F 0 {X) and (B) F t (X). In both cases A = 0.1 , T g = 2. 



field Fj (X) . In all cases the external modulation induces a periodic 
modulation of the mean field E with a period T a = 2jt/ w equal to 
the period of the modulation. At the same time, we have verified 
that, although the forcing term has zero average (i.e., it does not 
change the average input current), the average interspike interval 
is slightly self-adjusted and, what is more important, there is no 
evidence of locking between the modulation and the frequency of 
the single neurons. In other words, the behavior is similar to the 
spontaneous partial synchronization observed in van Vreeswijk 
(1996) (where the modulation is self-generated). 

Because of the unavoidable oscillations of the interspike inter- 
vals, it is necessary to identify the spike times with great care. In 
practice we integrate Equation (56) with a fixed time step Af, 
by employing a standard fourth-order Runge-Kutta integration 
scheme. At each time step we check if X 1 > 1, in which case we 
go one step back and adopt the Henon trick, which amounts to 
exchanging t and X 1 in the role of independent variable (Henon, 
1982). 

The linear stability analysis can be performed by linearizing 
the system (56), to obtain 

; dFQO) , 

x J = x + gd — A sin(<p)8(p, = 1, . . . , N, 

€ = p — at, 

p = -ap, 
Sip = 0; 

and by thereby estimating the corresponding Lyapunov spectrum. 

In the case of Fq and F\, we have always found that the 
Lyapunov spectrum scales as 1 /N 2 as theoretically predicted in 
the absence of external modulation (see Figure 5 for one instance 
of each of the two velocity fields). 

A similar agreement is also found for F3, where the Lyapunov 
spectrum scales as 1 /N 4 , exactly as in the absence of external forc- 
ing (see Figure 6). Analogous results have been obtained for the 
other velocity fields (data not shown), which confirm that the 
validity of the previous analysis extends to more complex dynam- 
ical regimes, as long as the membrane potentials are smoothly 
distributed. 

7. SUMMARY AND OPEN PROBLEMS 

In this paper we have discussed the linear stability of both fully 
synchronized and splay states in pulse-coupled networks of iden- 
tical oscillators. By following Abbott and van Vreeswijk (1993), 
we have obtained analytic expressions for the long-wavelength 
components of the Floquet spectra of the splay state for generic 
velocity fields and post synaptic potential profiles. The structure 
of the spectra depends on the smoothness of both the velocity 
field and the transmitted pulses. The smoother they are and the 
faster the eigenvalues decrease with the wavelength of the corre- 
sponding eigenvectors. In practice, while splay states arising in LIF 
neurons with 8-pulses have a finite degree of (in)stability along all 
directions, those emerging in analytic velocity fields have many 
exponentially small eigenvalues. These results have been derived 
in the mean field framework, where the system is assumed to be 
infinite. Although realistic neural networks are finite, the present 




k/N 



FIGURE 6 I Lyapunov spectra for neurons forced by an external 
periodic signal, we observe the scaling 1//V 4 for the continuous 
velocity field F 3 (X). The data refer to A = 0.1, T s = 2. 



analysis predicts correctly, even for finite systems, the stability of 
the eigenmodes associated to the fastest scales and the order of 
magnitude of the eigenvalues corresponding to slower time scales. 
Interestingly, the scaling behavior of the eigenvalues carries over 
to that of the Lyapunov exponents, when the network is periodi- 
cally forced, suggesting that our results have a relevance that goes 
beyond the highly symmetric solutions studied in this paper. 

Finally, we derived an analytic expression for the Floquet spec- 
tra for the fully synchronous state. In this case the exponents 
associated to the dynamics of the membrane potentials are all 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



February 2014 | Volume 8 | Article 8 | 10 



Olmi et al. 



Linear stability in networks of pulse-coupled neurons 



identical, as it happens for the diffusive coupling, but here the 
result is less trivial, due to the fact that one must take into account 
that arbitrarily close to the solution, the ordering of the neu- 
rons may be different. Moreover, the value of the (degenerate) 
Floquet exponent coincides with the evaporation exponent (van 
Vreeswijk, 1996; Pikovsky et al, 2001) whenever the pulses are 
sufficiently smooth, while for discontinuous pulses (like exponen- 
tial and 8-spikes) the equivalence is lost (see also di Volo et al., 
2013). 

For discontinuous velocity fields, another important property 
that has been confirmed by our analysis is the role of the ratio 
_R = N/(Tqu.) between the width of the single pulse (1/a) and the 
average interspike interval of the whole network (T = Tq/N). In 
fact, it turns out that the asynchronous regimes can be strongly 
stable along all directions only when R remains finite in the 
thermodynamic limit (and is possibly small). This includes the 
idealized case of 8 -like pulses, but also setups where the single 
pulses are so short that they can be resolved by the single neurons. 
Mathematically speaking, this result implies that the thermody- 
namic limit does not commute with the limit of a zero pulse- 
width. It would be interesting to check to what extent this prop- 
erty extends to more realistic models. A first confirmation result 
is contained in Pazo and Montbrio (2013), where the authors find 
a similar property in a network of Winfree oscillators. 

Among possible extensions of our analysis, one should defi- 
nitely mention the inclusion of delay in the pulse transmission. 
This generalization is far from trivial as it modifies the phase dia- 
gram of the possible states (see Bar et al., 2012 for a recent brief 
overview of the possible scenarios) and it complicates noticeably 
the stability analysis of the synchronized phase. An analytic treat- 
ment of this latter case is reported in Timme et al. (2002) for 
generic velocity fields and excitatory 8-pulses. 
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APPENDICES 

A. FOURIER COMPONENTS OF THE PHASE RESPONSE CURVE 

In this appendix we briefly outline the way the explicit expression 
of A n and B n , defined in Equation (23), can be derived in the large 
n limit for a velocity field F(X) that is either discontinuous, or 
continuous with discontinuous first derivatives at the border of 
the definition interval. 

The integration interval [0, 1] appearing in Equation (23) is 
splitted in n sub-intervals of length 1 /«, and the original equation 
can be rewritten as 



(A„ + iB„) 



n . 



k/n JClitny 

dy . 

\j(k-i)/ n G(y) 



(Al) 



For n sufficiently large we can assume that the variation of 1/ G(y) 
is quite limited within each sub-interval, and we can approximate 
the function as follows, up to the second order 



G(y) g + T Q F{y 0 ) 



T 0 F'(y 0 ) 
g + T 0 F(y 0 ) 



(y - yo) 



+ 



T 0 F'(y 0 ) 
g + T 0 F(y 0 ) 



TpF (y 0 ) 
2(£+T 0 F(y 0 )) 



(y-yoY 



where yp = (k — l)/n is the lower extremum of the nth sub- 
interval. 

By inserting these expansions into Equation (Al) and by 
performing the integration over the n sub-intervals, we can deter- 
mine an approximate expression for A„ and B n . The estimation of 
A n involves integrals containing cos(2jtny); it is easy to show that 
the integral over each sub-interval is zero if the integrand, which 
multiplies the cosinus term, is constant or linear in y; therefore 
the only non-zero terms are, 



pk/n 

J <k- r 



>(k-t)/n 
This allows to rewrite 



dy cos(2jt ny)y 
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2ti 2 h 3 



(A2) 
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where 



H 2 (x) 



TpF (x) 



(TpF'(x)) 1 

(£+r 0 F(x)) 3 2(g+T 0 F{x) 2 ) 



(A4) 



It is easy to verify that H2(x) admits an exact primitive and there- 
fore to perform the integral appearing in Equation (A3) and to 
arrive at the expression reported in Equation (28). 

The estimation of B n is more delicate, since now integrals con- 
taining sin(2n«y) are involved. The only vanishing integrals over 



the sub-intervals are those with a constant integrand multiplied 
by the sinus term and therefore the estimation of B„ reduces to 



'k — 1\ r k ' n 

!(k- \)/n 
*k/n 



dy sin(2jtny)y 
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and the non-zero integrals are 
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This allows to rewrite B„ as 
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We can then return to a continuous variable by rewriting (A8), up 
totheO(l/« 3 ),as 



1 

2tcm 



f 

.Jo 



Hi(x)dx + 



Hi(l)-Hi(0) 
In 



1 f l 

■- — - / H 2 (x)dx. 
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The expression Equation (29) is finally obtained by noticing that 
the primitive of H 2 (x) is Hi (x)/2, and that 



Jo 



H\(x)dx - 



1 



1 



(g+T o F(0)) (g+T 0 F(l)) 



For continuous velocity fields, B n = 0 so that, we can derive 
from Equation (26) an exact expression for the real part of the 
Floquet spectrum in the case of even L (for odd L the equivalent 
expression is given by Equation (31)) 

Wb ^o!^^, (A10) 



(2jt«)< 1 + 2 ) 
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A rigorous validation of the above formula would require going 
one order beyond in the l/n expansion of B n , a task that is 
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utterly complicated. In the specific case of the Quadratic Integrate 
and Fire neuron (or ©-neuron) F(X) = a — X(X — 1), it can be, 
however, analytically verified that B„ is exactly zero. Moreover, 
Equation (A10) is in very good agreement with the numeri- 
cally estimated Floquet spectra for two other continuous veloc- 
ity fields, namely F±(X) and -F2PO as shown in Figures 1, 3, 
respectively. As a consequence, it is reasonable to conjecture that 
Equation (29) is correct up to order C(l/n 4 ). 

B. EVAPORATION EXPONENT FOR THE LIF MODEL 

In this appendix we determine the (left and right) evaporation 
exponent for a synchronous state of a network of LIF neurons. 
This is done by estimating how the potential of a probe neu- 
ron, forced by the mean field generated by the network activity, 
converges toward the synchronized state. The stability analysis 
is performed by following the evolution of a perturbed probe 
neuron. Let us first consider an initial condition, where the syn- 
chronized cluster has just reached the threshold (X c = 1), while 
the probe neuron is lagging behind at a distance 8, . Such a distance 
is equivalent to a delay tn 



td 



Si 



F+(iy 



(All) 



where the subscript "+" means that the velocity field is estimated 
just after the pulses have been emitted. Over the time td, the 
potential of the cluster increases from the reset value 0 to 



hc=F+md= g® h , 



(A12) 



From now on (in LIF neurons), the distance decreases exponen- 
tially, reaching the value 



8/ = 8 c e T , 



(A13) 



after a period T. As a result, 



8/_P+(0) e _ r _ a + gE+ 



8, F+(l) a-l+gE+' 



(A14) 



The logarithm of the expansion factor gives the left evaporation 
exponent 

Let us now consider a probe neuron which precedes the syn- 
chronized cluster by an amount 8,. After a time T the distance 
becomes 



8 C = 8,e 



(A16) 



since no reset event has meanwhile occurred. Such a distance 
corresponds to a delay 



td' 



F~(l) 



(A17) 



where the subscript "— " means that the velocity has now to be 
estimated just before the pulse emission. By proceeding as before 
one obtains, 



8f = F-(0) _ r 
8, F-(l) 6 ■ 

so that the right evaporation exponent writes 



Al = ln 



a + gE- \ 
-1+gE-J 



T. 



(A18) 



(A19) 



It is easy to see that the left and right exponents differ if and only 
if E~ £+, i.e., if the pulses themselves are not continuous: this 
is, for instance, the case of exponential and 8 pulses. 
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